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I. INTRODUCTION 

The repulsive Hubbard model has long served as a sim- 
ple Hamiltonian to describe itinerant magnetism. Simi- 
larly, the attractive Hubbard model has been used to ex- 
plore qualitative features of the superconducting phase 
transition. [Q This model does not provide a microscopic 
model of the origin of pairing. Rather it is assumed that 
some other degrees of freedom, for example an electron- 
phonon coupling, have already provided the necessary at- 
traction. Recent quantum simulation studies have 
explored a number of features of the attractive Hubbard 
Hamiltonian, including a determination of the supercon- 
ducting transition temperature and its dependence on 
electron density, @] the detailed spatial structure of the 
pairing correlations, j^] the coexistence of a Bose-like spin 
ap with an otherwise degenerate Fermi gas of electrons, 
and deviations from Fermi-liquid behavior. Jlcfl Much 
analytic work has also been done, as reviewed in 

While the interplay of superconductivity and disorder 
has of course been extensively studied theoretically, fTl] ] 
considerably less is known numerically. The competition 
between the dephasing effect of impurity scattering and 
the dramatic manifestation of phase coherence in the zero 
resistance state gives rise to a set of challenging qual- 
itative questions. A quantitative understanding is also 
essential in order to model experiments like those which 
address the question of the possibility of a universal re- 
sistance in disordered superconducting films. |l2| 

Much recent theoretical |p^| , |l4T| and numerical 15 -li]] 
work on these issues has been done within the context 
of the "boson-Hubbard" model; that is, under the as- 
sumption that preformed Cooper pairs exist even in the 
non-superconducting state and the transition is driven 
by the loss of phase coherence, rather than the destruc- 
tion of the magnitude of the superconducting gap. This 
bosonic model should be the limit of the attractive Hub- 
bard Hamiltonian as the on-site interaction U becomes 
large. Despite the greater computational simplicity of 
the boson models, it is the fermion system which is of 
fundamental interest. 



In this paper we will explore the effect of random site 
energies in the attractive Hubbard Hamiltonian: 
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Here c\ a (c\ ) are operators which destroy (create) elec- 
trons of spin a on site i, so the first term in H describes 
the hopping of electrons between nearest neighbor sites 
on our 2d square lattice. \U\ is the on-site attraction, 
while \i and e\ are the chemical potential and random 
site energies, respectively. The e\ are chosen uniformly 
on [— V, +V]. In this paper, we will work exclusively at 
half-filling, {n^+riii) = 1. A related study of the current- 
current correlations and the behavior of the resistivity as 
a function of disorder strength and temperature away 
from half- filling is contained in |20|. 

In the absence of disorder, V = 0, considerable insight 
can be gained by considering the effect of a particle-hole 
transformation on the down electron operators, Cij «-» 
Citi Qi c^(-l)' ia!+ V The phase factor (— 
changes sign as one goes between the two sublattices of 
our (bipartite) square lattice. Under this transforma- 
tion the kinetic energy is invariant, while the interac- 
tion changes sign, \U\ —\U\. The chemical potential 
now couples to the z component of spin on each site, in- 
stead of to charge. In the absence of a chemical potential 
term, we have an exact mapping between the attractive 
and repulsive Hubbard Hamiltonians. Pair correlations 
in the attractive model map onto spin correlations in the 
xy plane of the repulsive model, while cdw correlations 
are associated with spin correlations in the z direction. 
Since the long range spin order in the ground state of the 
2D repulsive Hubbard model is rotationally invariant, we 
immediately conclude that pairing and cdw correlations 
coexist in the ground state of the half-filled attractive 
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model. This mapping can also be used to discuss the 
effect of doping, as will be described in section 5. 

In this manuscript we consider non-zero disorder, V =/= 
0. Our basic conclusions are: (i) The addition of site 
disorder breaks the degeneracy between superconduct- 
ing and charge-density wave states. Equal time density- 
density correlations exhibit a rapid suppression of their 
staggered pattern, while pair-pair correlations remain ro- 
bust, (ii) Our data is consistent with an immediate de- 
struction of long range charge order, while supercon- 
ducting correlations appear to persist out to a finite 
V c » 1.5. (iii) The behavior of superconducting corre- 
lations at \U\ — 4 are in approximate agreement with 
that of the appropriate magnetic correlations in a strong 
coupling model. 

An outline of the remainder of this paper is as follows: 
In section 2 we give a brief description of the simulation 
and of the observables used to characterize the ground 
state correlations. In section 3 we present results for 
local pair and charge correlations. Section 4 provides 
a finite size scaling analysis of this data to determine 
the existence of long range order. We also discuss in 
some detail the distribution of measurements for differ- 
ent disorder realizations. Section 5 discusses results for 
the strong coupling version of this model, the antiferro- 
magnetic Heisenberg model in a random magnetic field. 
A summary is presented in section 6. 



II. BRIEF DESCRIPTION OF THE SIMULATION 

We will use the "determinant" monte carlo algorithm 
pl[ for our numerical work. In this approach the par- 
tition function is written down as a path integral by 
discretizing the inverse temperature (3 and using the 
Trotter approximation p^psj ] to break up the exponen- 
tial of the kinetic and potential pieces of the Hamilto- 
nian. The interaction term is decoupled with a discrete 
Hubbard- Stratonovich transformation p4| . The result- 
ing trace over the fermion operators is over the exponen- 
tial of quadratic forms, and so can be done analytically. 
The result is an expression for the partition function 
which is the sum over all values of the discrete Hubbard- 
Stratonovich field of a summand which is the product 
of two determinants, one each arising from the spin up 
and the spin down degrees of freedom. Because both 
the random site energies and the auxiliary field couple to 
the charge, the two determinants are identical and hence 
their product is positive. There is no fermion sign prob- 
lem in these simulations, even away from half filling. 

The matrices M a whose determinants give the weight 
of a particular auxiliary field configuration are simply 
related to the equal time fermion Green function: Gy = 

(ciCj) = M^j 1 . Observables are measured by expressing 
them (using Wick's theorem) in terms of the appropriate 
sums and products of G. Of particular interest to us here 
are the equal time charge and pair correlations, 
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and their associated structure factors, 
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FIG. 1. Two possible real space configurations of the elec- 
trons on the lattice are shown. The pair structure factor at- 
tains its largest value N/4 for both, while the cdw structure 
factor is large only for (a). 

It is useful to look at a few strong coupling snapshots 
of possible real space electron distributions to gain a pre- 
liminary insight into these correlation functions and the 
effect of site disorder. Consider the strong coupling limit 
when all electrons in the system are paired. A typical low 
energy state consists then of a lattice with each site either 
empty or doubly occupied. A configuration in which the 
doubly occupied sites alternate with empty sites (Fig. la) 
has a lower energy than one in which doubly occupied or 
empty sites are adjacent by an amount AE cx —t 2 /U to 
second order in the hopping t. (In the language of the 
repulsive model onto which the attractive model maps 
via the particle-hole transformation discussed above, this 
energy lowering stabilizes antifcrromagnetism over ferro- 
magnetism at half-filling, and is proportional to the ex- 
change constant J.) The charge density and pair struc- 
ture factors defined in Eq. 3 take on their maximal values 
(N/4 for Vir and N for S cdw ). 

Now consider the effect of site disorder. It will not 
lead to the breaking of the pairs, but, when the site en- 
ergies exceed cx t 2 /U, it will change the sites on which 
the pairs prefer to reside (Fig. lb). Note that the pair 
structure factor is still large for such a disordered config- 
uration of pairs, since contributions to it depend only on 
finding doubly occupied sites and empty sites somewhere 
in the lattice to which to hop. However, the phases in 
the charge structure factor make it extremely sensitive 
to the destruction of the original staggered pattern. This 
rough argument suggests that pair order will be more ro- 
bust to randomness in the site energy than will charge 
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ordering. Of course, on general grounds we also expect 
a term in the Hamiltonian which couples directly to the 
charge to have the greatest effect on the associated charge 
correlations. Indeed, as discussed by Anderson [^6) non- 
magnetic impuritites are not expected to destroy super- 
conductivity, since one can still pair appropriate eigen- 
states of the single particle Hamiltonian which includes 
the randomness. Even when the disorder is large enough 
to localize these eigenstates, it has been suggested that 
superconductivity survives. (27) 

This overview captures the essence of how disorder af- 
fects our system. In the next section, we will make this 
qualitative picture more precise. 



III. LOCAL CORRELATIONS 
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We begin by showing some results for the disorder de- 
pendence of local quantities. The kinetic energy (k) = 



-jf L<i,j) CT ( c L c jT + c J v Ci ' J ) is shown in Fi S- 2 - As we sna11 
see below, superconductivity vanishes around V c w 1.5. 
The kinetic energy shows no special signal at this tran- 
sition. Of course a measure of local electron hopping 
like (k) does not have to vanish at an insulating phase 
transition, pal 




FIG. 2. The kinetic energy as a function of disorder 
strength. 



FIG. 3. The double occupancy rate as a function of disor- 
der strength. 
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FIG. 4. Charge correlations as a function of site separation 
for different disorder strengths. Successive disorder strengths 
have been offset vertically for clarity. The horizontal lines in- 
dicate the c(j - 1)=0 axis for each successive disorder strength. 
(The numbers labeling the vertical axis correspond to the 
V=0.00 case.) 
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FIG. 5. Pair correlations as a function of site separation 
for different disorder strengths. 



a function of lattice separation j — 1 for different disor- 
der strengths. (Successive disorder strengths have been 
offset vertically for clarity.) The lattice size is 8x8, in- 
verse temperature (3 = 10 and \U\ = 4. The oscillatory 
character of the charge correlations is indicative of cdw 
ordering. At V = these correlations extend over the en- 
tire lattice, that is, the correlation length £ exceeds the 
linear lattice dimension. However as V is turned on the 
correlations go to zero. 

Fig. 5 shows the analogous plot for the pair correla- 
tion function p s (j — 1) . The pair correlations remain un- 
changed for weak V , then are eventually suppressed for 
sufficiently large randomness. We see the robustness of 
the pair correlations as compared to charge-charge cor- 
relations in this plot, which has a much greater range of 
disorder strengths, V, than does Fig. 4. 

The Fourier transforms of these real space correlation 
functions are shown as a function of disorder in Fig. 6. 
The degeneracy between charge and pair correlations is 
evident in the absence of randomness, V = 0. As was 
seen in Figs. 4 and 5, nonzero site disorder more rapidly 
destroys the charge density wave than the pair correla- 
tions. 
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FIG. 6. Charge and pair structure factors as a function of 
disorder strength. 



IV. THE DESTRUCTION OF LONG RANGE 
ORDER 

To determine whether ground-state long range order 
exists in our system we need to do a finite size scal- 
ing analysis. As has been discussed ||3l]]34 within the 
context of the repulsive Hubbard model, spin-wave the- 
ory (3^] predicts that on a 2D lattice of size N=LxL the 
charge structure factor and correlation function at largest 
separation should behave as 
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Scdw = m /3 + a/L 



c(L/2,L/2) = m 2 /3 + b/L 



(4) 



in the ordered phase. Similar results are valid for the 
pair correlations. Here m is the order parameter. Thus, 
in the ordered phase a plot of the scaled structure factor 
versus 1/L should be a straight line with a nonzero in- 
tercept giving m 2 /3. We will always choose the inverse 
temperature /3 sufficiently large that we are effectively at 
T = for our finite lattices. 



In the repulsive Hubbard model, random site energies 
have a fundamental qualitative effect on the double occu- 
pancy rate (jiji-fijj.) since they compete with the repulsive 
interaction and disturb local moment formation. In the 
attractive model, we expect site disorder to have a much 
less dramatic effect and indeed that is seen to be the case 
in Fig. 3. 

Longer range charge-charge correlations are dramati- 
cally affected by site disorder. Fig. 4 shows c(j — 1) as 



4 



0.2 



0.1 



I I I I 

. (a) V = 0.00 










_ (b) V = 0.25 


n 






. (c) V = 0.50 









_ (d) V=0.75 














- 




A 







0.1 



0.2 



1/L 



FIG. 7. Finite size scaling plots for the charge correlations: 
(a) V=0.0; (b) V=0.25; (c) V=0.50; and (d) V=0.75. The 
straight lines are least-squares fits to the data. Error bars 
(not shown) on the V=0.25 case are consistent with a zero 
intercept, while error bars on the V=0.00 case are not con- 
sistent with a zero intercept. The V=0.00 (clean case) error 
bars are much smaller than in the disordered V=0.25 case, 
since there is no disorder averaging required in the clean case. 
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FIG. 9. Histogram of values of the (a) interaction and (b) 
kinetic energies for different disorder realizations. Each plot 
represents data from one hundred disorder realizations. The 
average values are shown by the vertical lines. 
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FIG. 8. Finite size scaling plots for the pair correlations: 
(a) V=0.0; (b) V=0.50; (c) V=1.00; and (d) V=1.50. The 
straight lines are least-squares fits to the data. 



We will conclude this section with a discussion of dis- 
order averaging, since while the other aspects of our sim- 
ulation are identical to those long reported for determi- 
nant monte carlo, little is known about what happens 
when randomness is included. In Fig. 9 we show some 
histograms of the interaction and kinetic energies for an 
8x8 spatial lattice at inverse temperature /3 = 10. We 
have chosen \U\ = 4 and V = 1. We see that these quan- 
tities have a fairly sharp distribution, that is, the energy 
is not too sensitive to the detailed disorder realization; 
the widths of the distributions are less than 5 % of the av- 
erage values. The error bars associated with realization 
to realization fluctuations in these quantities are roughly 
ten times the statistical uncertainties in a run consisting 
of 1000 warm-up sweeps and 5000 measurement sweeps 
for a single realization. 



Fig. 7 shows the result of this analysis for the charge 
correlations. The interaction strength is \U\ = 4. Only in 
the clean system at V = is a nonzero order parameter 
m obtained. However, as seen in Fig. 8 the pair field 
order parameter remains nonzero out to approximately 
V = V C * 1.5. |§ 
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FIG. 10. Histogram of values of the (a) charge and (b) pair 
structure factors for different disorder realizations. Each plot 
represents data from 100 disorder realizations. The average 
values are shown by the vertical lines. The charge structure 
factor, S c dw, has an apparently non-gaussian form. 

On the other hand, the fluctuations in quantities which 
measure long range correlations are, as expected, much 
larger. In Fig. 10 we show histograms of the charge and 
pair structure factors for an 8x8 spatial lattice at inverse 
temperature (3 = 10, with \U\ = 4 and V = 1. We see 
that the widths of the peaks are of the same order of mag- 
nitude as the average values of the respective structure 
factors. 

We note that practical difficulties limit the number of 
disorder realizations, and consequently the reduction in 
fluctuations due to disorder averaging, that one can at- 
tain. A single disorder realization for an 8x8 system at 
P = 10 takes over one-thousand cpu minutes on a fast 
workstation. This computational difficulty of these QMC 
calculations precludes disorder averaging over thousands 
of realizations, as has been done in the spin glass litera- 
ture. The non-gaussian nature of the distributions (e.g., 
Fig. 10a) of course raises difficult questions about how 
to do the averaging and how to estimate error bars cor- 
rectly. However, if one goes ahead and employs the usual 
methods of getting error bars based on an assumption of 
a gaussian distribution, then averaging over 20-100 dis- 
order realizations reduces the statistical errors to about 
the same level as the statistical errors associated with the 
monte carlo sampling. This is what we have done in the 
data reported in this paper. 



V. THE RANDOM FIELD HEISENBERG MODEL 

As we have discussed, the attractive (— U) Hubbard 
model can be mapped onto the repulsive model, which in 
turn at strong coupling can be mapped onto a quantum 
spin 1/2 antiferromagnetic Heisenberg Hamiltonian. In 
the absence of disorder, the behavior of the associated 
classical spin models has exhibited considerable analo- 
gies to the original —U Hubbard model. |2| Here we de- 
sire to see if similar connections can usefully be made 
between the disordered, —U Hubbard model and the as- 
sociated classical model - the random field Heisenberg 
model. However, the problem of the random field Heisen- 
berg model is an extremely difficult one in it own right. 
We emphasize that we are attempting only qualitative 
contact with the attractive Hubbard model simulations 
here. 

We begin by reviewing the results in a uniform mag- 
netic field, since the comparison will be useful in dis- 
cussing the case of a random field. Similar results were 
presented in However, here we present some addi- 
tional plots which help more precisely characterize the 
nature of the ordered phase. In the absence of a field, 
the continuous symmetry of the model assures us that 
in 2D there can be no true long range order except in 
the ground state T — 0. Q If a field h z — is ap- 
plied, the spins tend to lie down in the xy plane, because 
then they can tilt upwards in the z direction and take 
advantage of the field energy without costing as much 
exchange energy J as if they were antiferromagnetically 
aligned in the z direction. Thus, the antiferromagnetic 
Heisenberg model in a uniform magnetic field is argued 
to be in the universality class of the XY model, with a 
finite temperature Kosterlitz-Thouless phase transition 
into a state with power law decay of the correlation func- 
tions. In the language of the attractive model, doping 
breaks the cdw-pair degeneracy, and off half-filling one 
has a finite temperature phase transition into a purely 
superconducting state. (t]] 

Let us define the antiferromagnetic structure factors in 
the different spin directions a — x, y, z as the appropriate 
sums of the correlations of the a component of spin on 
different sites, 

s aa = i^( Sa aK(i))(-i) |j - 11 . (5) 
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FIG. 11. Results for the structure factors as a function of 
temperature for a 20x20 lattice, (a) S xx + S yy ; (b) 2S ZZ . The 
line at T=0.67 is at the temperature used in the finite size 
scaling plot (Fig. 12) below. 
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FIG. 12. Results for the xy structure factors as a function 
of uniform field h z for different lattice sizes. 

Simulation results for the classical antiferromagnetic 
Heisenberg model in a uniform magnetic field are shown 
in Figs. 11-13. The temperature dependence is shown 
first in Fig. 11 on a fixed lattice size. To set a possible 
scale of T, note that T KT = 0.725 for the XY model. 
p5| This is consistent with the temperature at which 
the structure factor swings upward in Fig. 11. We note 
that a small uniform field, h z — 1.0, enhances S xx and 
Syy substantially, and dramatically reduces S zz . Fig. 12 
shows the field dependence of the xy structure factors 
at a fixed T = 0.67 for different lattices. @ When h z 
is nonzero, there is a significant size dependence of the 
structure factor even at nonzero temperature, which sug- 
gests that the presence of a field may indeed make the 
system order at finite temperature. 

Of course this data is only suggestive. A careful fi- 
nite size scaling analysis would be needed to pin down 
whether T c = or T c ^ 0. To illustrate one of the is- 
sues involved, we note that even if a phase transition 
occurs only at T = 0, structure factors will show sig- 
nificant size dependence up to the lattice size, L rs £. 
In the Heisenberg model, the corrrelation length £ = 
Qexp(27rJ/r)/(l + 2vr J/T) where C ( « 0.01. @ Here, 
with T = 0.67, the correlation length £ m 11. The growth 
in the structure factors S aa with lattice size at h z = is 
consistent with this value of £, and hence with T c = 0, as 
should be the case. The increased growth of S aa with lat- 
tice size at nonzero h z could merely reflect a larger £ (but 
T c still zero) or a finite T c . The structure factors, S aa , 
continue to show the effect of this finite £ throughout the 
range of lattice sizes used in the simulations illustrated in 
Fig. 12. Hence, the value of £ in this model is consistent 
with the results shown in Fig. 12 being at a temperature 
where correlations have started to form across the lattice. 
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32x32 Lattice 




FIG. 13. The averages of the squares of different spin com- 
ponents C aa = {s'a} as a function of uniform field h z . For 
h z = 0, C X x = Cyy — C Z z = 1/3 from rotational invariance. 
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We note that while the picture of the spins "lying" 
down in the xy plane is a physically appealing one, it 
is not in fact so accurate a view quantitatively. Fig. 13 
shows the square of the spin components as a function 
of field strength. We show only one lattice size, 32x32, 
since there is no dependence on lattice size for such a 
local quantity. It is apparent that even for h z = 2.0 
where Fig. 12 shows a very significant enhancement of 
the xy plane structure factor, the spins locally still point 
almost as much in the z direction as for h z = 0.0. Indeed 
for larger fields, the spins start to align ferromagncti- 
cally with the field and (s 2 z ) exceeds (s 2 ) — (Sy). How 
do we reconcile this with the results for the structure 
factor? Apparently, the uniform field h z has a large ef- 
fect on the long range spin correlations, and very little 
effect on the short range correlations. Thus while indi- 
vidual spins still rotate significantly in the z direction, 
the antiferromagnetic correlations in the z direction be- 
tween different spins are destroyed. (In fact they become 
ferromagnetic.) At the same time, long range antifer- 
romagnetic correlations in the plane are enhanced, even 
though the individual spins are not really "lying down" 
in the plane. 



o 
o 



O 
O 



o 
o 



h z = 0.0 o 
Ah =1 .0 □ 




(b) 
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FIG. 14. Results for the structure factors as a function of 
temperature for a 20x20 lattice in both a zero field (h z = 0.0) 
and random field (Ah z — 1.0). (a) S xx + S yy ; (b) 2S ZZ . The 
line at T=0.3 is at the temperature used in the finite size 
scaling plots (Fig. 15) below. 
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FIG. 15. Results for the structure factors at T=0.3 as a 
function of random field strength Ah z for different lattice 
sizes, (a) S xx + S yy ; (b) 2S ZZ . 

We now turn to the results of simulations for random 
magnetic fields h z (j), which is the situation to which 
our random site energy attractive Hubbard model cor- 
responds. Figs. 14-15 show the analogous plots to those 
of the uniform field case. Again for low T a random field 
enhances S xx and S yy , although the effect is much less 
dramatic than seen in Fig. 11a. The suppression of S zz 
is more substantial, but again less decisive than with a 



uniform field, Fig. lib. The random-field results shown 
in Fig. 15 were obtained at T=0.3 while the results in 
the corresponding uniform-field Fig. 12 were obtained at 
T=0.67, the temperature at which previously reported 
uniform-field Heisenberg work was done. D For interme- 
diate temperatures like T — 0.67 it was found that S zz 
exceeds, by a small amount, S xx and S yy for small dis- 
order. It is not clear to us why this occurs. As seen 
by comparing Figs, lib and 14b, S zz behaves quite dif- 
ferently when the system is subjected to a uniform as 
opposed to a random field along the z axis. A uniform 
field suppresses S zz much more strongly than a random 
field. Nevertheless, in both cases the z axis structure fac- 
tor is suppressed with respect to the xy plane structure 
factors at low enough temperatures. 

In Fig. 15, at T=0.3, we see that increasing the disor- 
der suppresses the structure factor in the z direction, the 
analogue of the cdw structure factor (Fig. 15b), while the 
structure factors in the xy plane, the analogue of the pair 
structure factor, remain robust to randomness (Fig. 15a). 
Hence, we see the same general results as seen in the 
disordered, attractive Hubbard model. The analogy is 
somewhat limited, since we do not see the destruction 
of the xy plane structure factor (Fig. 15a) as disorder is 
increased, the analogue of the destruction of supercon- 
ductivity (Fig. 6). Determining the critical properties of 
the random field Heisenberg model is a problem far be- 
yond the scope of this paper. However, we note that a 
finite V c for S xx and S yy is not forbidden, for example, 
by an Imry-Ma argument, since although the order pa- 
rameter in this case has continuous symmetr y, i t is not 
conjugate to the random disordered h z field. ]3q| 



VI. CONCLUSIONS 

We have studied the effect of random site energies on 
charge density wave and superconducting correlations in 
the 2D attractive Hubbard Hamiltonian at half-filling. 
We find an immediate suppression of cdw order with a 
random one-body potential V , whereas the pair order is 
relatively robust. This is consistent with a strong cou- 
pling picture of the real space configurations of low en- 
ergy states of pairs as disorder is turned on. 

After showing the effect of V on the spatial correla- 
tions, we performed a finite size scaling analysis of the 
structure factors and determined that long range cdw is 
destroyed for V c w 0.0 while superconducting order has 
V c « 1.5 which is roughly the energy scale 4t 2 /U rj 1 
which stabilizes pairing. It is interesting that V c is con- 
siderably less than that required for the destruction of 
superconductivity in the doped system. p(jf 

We discussed briefly the distribution of different mea- 
surements as the disorder configuration was modified. 
This enabled us to determine a reasonable number of 
disorder realizations to average over. Finally, the strong 
coupling picture of the classical Heisenberg Hamiltonian 
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in a random external field h z (j) was shown to exhibit 
some similarities with the disordered attractive Hubbard 
model results at \U\ — 4.0. 

There are a number of further questions we wish to ex- 
plore. First, we always performed our simulations at an 
inverse temperature f3 sufficiently large so that we were 
in the ground state of the finite size lattices we were 
studying. The temperature dependence of various quan- 
tities would be interesting to determine. We would also 
like to look more at the non-equal time, dynamical re- 
sponse of the system. In particular, the issue of how the 
superconducting gap fills in with disorder is an interest- 
ing one. This work requires analytic continuation of the 
imaginary time Green's function, a task which should be 
challenging, especially in the presence of disorder. 
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